#include "antsUtilities.h"
#include "antsAllocImage.h"
#include <algorithm>
#include "ReadWriteData.h"

#include "itkAffineTransform.h"
#include "itkCSVArray2DDataObject.h"
#include "itkCSVArray2DFileReader.h"
#include "itkCSVNumericObjectFileWriter.h"
#include "itkImage.h"
#include "itkLabelGeometryImageFilter.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkResampleImageFilter.h"
#include "itkTransformFileWriter.h"

#include "itkLabelPerimeterEstimationCalculator.h"
#include "itkLabelMap.h"
#include "itkLabelImageToShapeLabelMapFilter.h"
#include "itkShapeLabelMapFilter.h"
#include "itkShapeLabelObject.h"

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <ostream>
#include <sstream>
#include <string>

namespace ants
{
template <unsigned int ImageDimension>
int
LabelGeometryMeasures(int argc, char * argv[])
{
  using LabelType = unsigned int;
  using LabelImageType = itk::Image<LabelType, ImageDimension>;

  using RealType = float;
  using RealImageType = itk::Image<RealType, ImageDimension>;

  typename LabelImageType::Pointer labelImage = LabelImageType::New();
  ReadImage<LabelImageType>(labelImage, argv[2]);

  typename RealImageType::Pointer intensityImage = RealImageType::New();
  bool                            intensityImageUsed = false;
  if (argc > 3)
  {
    std::string intensityImageFile = argv[3];
    ConvertToLowerCase(intensityImageFile);
    if ( intensityImageFile.compare("na") != 0 && intensityImageFile.compare("none") != 0 )
    {
      if (!ReadImage<RealImageType>(intensityImage, intensityImageFile.c_str()))
      {
        std::cout << "Input intensity image " << intensityImageFile << " could not be read." << std::endl;
        return EXIT_FAILURE;
      }
      intensityImageUsed = true;
    }
  }

  bool outputCSVFormat = false;
  if (argc > 4)
  {
    outputCSVFormat = true;
  }
  using FilterType = itk::LabelGeometryImageFilter<LabelImageType, RealImageType>;
  typename FilterType::Pointer filter = FilterType::New();
  filter->SetInput(labelImage);
  if (intensityImageUsed)
  {
    filter->SetIntensityInput(intensityImage);
  }
  filter->CalculatePixelIndicesOff();
  filter->CalculateOrientedBoundingBoxOff();
  filter->CalculateOrientedLabelRegionsOff();
  // These generate optional outputs.
  //   filter->CalculatePixelIndicesOn();
  //   filter->CalculateOrientedBoundingBoxOn();;
  //   filter->CalculateOrientedLabelRegionsOn();

  filter->Update();

  using LabelObjectType = itk::ShapeLabelObject<LabelType, ImageDimension>;
  using LabelMapType = itk::LabelMap<LabelObjectType>;

  // convert the image in a collection of objects
  using ConverterType = itk::LabelImageToShapeLabelMapFilter<LabelImageType, LabelMapType>;
  typename ConverterType::Pointer converter = ConverterType::New();
  converter->SetInput(labelImage);

  using ValuatorType = itk::ShapeLabelMapFilter<LabelMapType>;
  typename ValuatorType::Pointer valuator = ValuatorType::New();
  valuator->SetInput(converter->GetOutput());

  valuator->Update();

  typename LabelMapType::Pointer labelMap = valuator->GetOutput();

  //   typedef itk::LabelPerimeterEstimationCalculator<LabelImageType> AreaFilterType;
  //   typename AreaFilterType::Pointer areafilter = AreaFilterType::New();
  //   areafilter->SetImage( labelImage );
  //   areafilter->SetFullyConnected( false );
  //   areafilter->Compute();

  typename FilterType::LabelsType allLabels = filter->GetLabels();
  std::sort(allLabels.begin(), allLabels.end());

  if (outputCSVFormat)
  {
    typename FilterType::LabelsType::iterator allLabelsIt;

    std::vector<std::string> columnHeaders;

    columnHeaders.emplace_back("Label");
    columnHeaders.emplace_back("VolumeInVoxels");
    columnHeaders.emplace_back("SurfaceAreaInMillimetersSquared");
    columnHeaders.emplace_back("Eccentricity");
    columnHeaders.emplace_back("Elongation");
    columnHeaders.emplace_back("Orientation");
    columnHeaders.emplace_back("Centroid_x");
    columnHeaders.emplace_back("Centroid_y");
    if (ImageDimension == 3)
    {
      columnHeaders.emplace_back("Centroid_z");
    }
    columnHeaders.emplace_back("AxesLength_x");
    columnHeaders.emplace_back("AxesLength_y");
    if (ImageDimension == 3)
    {
      columnHeaders.emplace_back("AxesLength_z");
    }
    columnHeaders.emplace_back("BoundingBoxLower_x");
    columnHeaders.emplace_back("BoundingBoxUpper_x");
    columnHeaders.emplace_back("BoundingBoxLower_y");
    columnHeaders.emplace_back("BoundingBoxUpper_y");
    if (ImageDimension == 3)
    {
      columnHeaders.emplace_back("BoundingBoxLower_z");
      columnHeaders.emplace_back("BoundingBoxUpper_z");
    }

    if (filter->GetIntensityInput())
    {
      columnHeaders.emplace_back("IntegratedIntensity");
      columnHeaders.emplace_back("WeightedCentroid_x");
      columnHeaders.emplace_back("WeightedCentroid_y");
      if (ImageDimension == 3)
      {
        columnHeaders.emplace_back("WeightedCentroid_z");
      }
    }

    std::vector<std::string> rowHeaders;
    for (allLabelsIt = allLabels.begin(); allLabelsIt != allLabels.end(); allLabelsIt++)
    {
      if (*allLabelsIt == 0)
      {
        continue;
      }
      std::ostringstream convert; // stream used for the conversion
      convert << *allLabelsIt;    // insert the textual representation of 'Number' in the characters in the stream
      rowHeaders.push_back(convert.str()); // set 'Result' to the contents of the stream
    }

    vnl_matrix<double> measures(allLabels.size() - 1, columnHeaders.size() - 1);

    unsigned int rowIndex = 0;
    for (allLabelsIt = allLabels.begin(); allLabelsIt != allLabels.end(); allLabelsIt++)
    {
      if (*allLabelsIt == 0)
      {
        continue;
      }

      unsigned int columnIndex = 0;
      //      measures( rowIndex, columnIndex ) = static_cast< double >( *allLabelsIt );
      measures(rowIndex, columnIndex++) = filter->GetVolume(*allLabelsIt);

      const LabelObjectType * labelObject = labelMap->GetLabelObject(*allLabelsIt);

      measures(rowIndex, columnIndex++) = labelObject->GetPerimeter();
      measures(rowIndex, columnIndex++) = filter->GetEccentricity(*allLabelsIt);
      measures(rowIndex, columnIndex++) = filter->GetElongation(*allLabelsIt);
      measures(rowIndex, columnIndex++) = filter->GetOrientation(*allLabelsIt);
      measures(rowIndex, columnIndex++) = filter->GetCentroid(*allLabelsIt)[0];
      measures(rowIndex, columnIndex++) = filter->GetCentroid(*allLabelsIt)[1];
      if (ImageDimension == 3)
      {
        measures(rowIndex, columnIndex++) = filter->GetCentroid(*allLabelsIt)[2];
      }
      measures(rowIndex, columnIndex++) = filter->GetAxesLength(*allLabelsIt)[0];
      measures(rowIndex, columnIndex++) = filter->GetAxesLength(*allLabelsIt)[1];
      if (ImageDimension == 3)
      {
        measures(rowIndex, columnIndex++) = filter->GetAxesLength(*allLabelsIt)[2];
      }

      unsigned int arrayIndex = 0;

      measures(rowIndex, columnIndex++) = filter->GetBoundingBox(*allLabelsIt)[arrayIndex++];
      measures(rowIndex, columnIndex++) = filter->GetBoundingBox(*allLabelsIt)[arrayIndex++];
      if (ImageDimension == 3)
      {
        measures(rowIndex, columnIndex++) = filter->GetBoundingBox(*allLabelsIt)[arrayIndex++];
      }
      measures(rowIndex, columnIndex++) = filter->GetBoundingBox(*allLabelsIt)[arrayIndex++];
      measures(rowIndex, columnIndex++) = filter->GetBoundingBox(*allLabelsIt)[arrayIndex++];
      if (ImageDimension == 3)
      {
        measures(rowIndex, columnIndex++) = filter->GetBoundingBox(*allLabelsIt)[arrayIndex++];
      }

      if (filter->GetIntensityInput())
      {
        measures(rowIndex, columnIndex++) = filter->GetIntegratedIntensity(*allLabelsIt);
        measures(rowIndex, columnIndex++) = filter->GetWeightedCentroid(*allLabelsIt)[0];
        measures(rowIndex, columnIndex++) = filter->GetWeightedCentroid(*allLabelsIt)[1];
        if (ImageDimension == 3)
        {
          measures(rowIndex, columnIndex++) = filter->GetWeightedCentroid(*allLabelsIt)[2];
        }
      }
      rowIndex++;
    }

    using WriterType = itk::CSVNumericObjectFileWriter<double, 1, 1>;
    WriterType::Pointer writer = WriterType::New();
    writer->SetFileName(argv[4]);
    writer->SetColumnHeaders(columnHeaders);
    writer->SetRowHeaders(rowHeaders);
    writer->SetInput(&measures);
    try
    {
      writer->Write();
    }
    catch (const itk::ExceptionObject & exp)
    {
      std::cerr << "Exception caught!" << std::endl;
      std::cerr << exp << std::endl;
      return EXIT_FAILURE;
    }
  }
  else
  {
    typename FilterType::LabelsType::iterator allLabelsIt;
    //   std::cout << "Number of labels: " << labelGeometryFilter->GetNumberOfLabels() << std::endl;
    //   std::cout << "Label geometry measures." << std::endl;
    std::cout << std::left << std::setw(7) << "Label" << std::left << std::setw(10) << "Volume(voxels)" << std::left
              << std::setw(15) << "SurfArea(mm^2)" << std::left << std::setw(15) << "Eccentricity" << std::left
              << std::setw(15) << "Elongation" << std::left << std::setw(15) << "Orientation" << std::left
              << std::setw(30) << "Centroid" << std::left << std::setw(30) << "Axes Length" << std::left
              << std::setw(30) << "Bounding Box";
    if (filter->GetIntensityInput())
    {
      std::cout << std::left << std::setw(20) << "Integrated Int." << std::left << std::setw(30) << "Weighted Centroid";
    }
    std::cout << std::endl;
    for (allLabelsIt = allLabels.begin(); allLabelsIt != allLabels.end(); allLabelsIt++)
    {
      if (*allLabelsIt == 0)
      {
        continue;
      }
      std::cout << std::setw(7) << *allLabelsIt;
      std::cout << std::setw(10) << filter->GetVolume(*allLabelsIt);

      const LabelObjectType * labelObject = labelMap->GetLabelObject(*allLabelsIt);

      std::cout << std::setw(15) << labelObject->GetPerimeter();
      std::cout << std::setw(15) << filter->GetEccentricity(*allLabelsIt);
      std::cout << std::setw(15) << filter->GetElongation(*allLabelsIt);
      std::cout << std::setw(15) << filter->GetOrientation(*allLabelsIt);

      std::stringstream oss;
      oss << filter->GetCentroid(*allLabelsIt);
      std::cout << std::setw(30) << (oss.str()).c_str();
      oss.str("");

      oss << filter->GetAxesLength(*allLabelsIt);
      std::cout << std::setw(30) << (oss.str()).c_str();
      oss.str("");

      oss << filter->GetBoundingBox(*allLabelsIt);
      std::cout << std::setw(30) << (oss.str()).c_str();
      oss.str("");

      //     std::cout << filter->GetMajorAxisLength( *allLabelsIt ) << "\t";
      //     std::cout << filter->GetMinorAxisLength( *allLabelsIt ) << "\t";
      if (filter->GetIntensityInput())
      {
        oss << filter->GetIntegratedIntensity(*allLabelsIt);
        std::cout << std::setw(20) << (oss.str()).c_str();
        oss.str("");

        oss << filter->GetWeightedCentroid(*allLabelsIt);
        std::cout << std::setw(30) << (oss.str()).c_str();
        oss.str("");
      }
      std::cout << std::endl;
    }
  }

  return EXIT_SUCCESS;
}

// entry point for the library; parameter 'args' is equivalent to 'argv' in (argc,argv) of commandline parameters to
// 'main()'
int
LabelGeometryMeasures(std::vector<std::string> args, std::ostream * itkNotUsed(out_stream))
{
  // put the arguments coming in as 'args' into standard (argc,argv) format;
  // 'args' doesn't have the command name as first, argument, so add it manually;
  // 'args' may have adjacent arguments concatenated into one argument,
  // which the parser should handle
  args.insert(args.begin(), "LabelGeometryMeasures");

  int     argc = args.size();
  char ** argv = new char *[args.size() + 1];
  for (unsigned int i = 0; i < args.size(); ++i)
  {
    // allocate space for the string plus a null character
    argv[i] = new char[args[i].length() + 1];
    std::strncpy(argv[i], args[i].c_str(), args[i].length());
    // place the null character in the end
    argv[i][args[i].length()] = '\0';
  }
  argv[argc] = nullptr;
  // class to automatically cleanup argv upon destruction
  class Cleanup_argv
  {
  public:
    Cleanup_argv(char ** argv_, int argc_plus_one_)
      : argv(argv_)
      , argc_plus_one(argc_plus_one_)
    {}

    ~Cleanup_argv()
    {
      for (unsigned int i = 0; i < argc_plus_one; ++i)
      {
        delete[] argv[i];
      }
      delete[] argv;
    }

  private:
    char **      argv;
    unsigned int argc_plus_one;
  };
  Cleanup_argv cleanup_argv(argv, argc + 1);

  // antscout->set_stream( out_stream );

  if (argc < 3)
  {
    std::cout << "Usage 1: " << argv[0] << " imageDimension labelImage [intensityImage] [csvFile]" << std::endl;
    std::cout << "To output to csvFile without an intensity image, use \"none\" or \"na\" as the third arg" << std::endl;
    if (argc >= 2 && (std::string(argv[1]) == std::string("--help") || std::string(argv[1]) == std::string("-h")))
    {
      return EXIT_SUCCESS;
    }
    return EXIT_FAILURE;
  }

  switch (std::stoi(argv[1]))
  {
    case 2:
    {
      LabelGeometryMeasures<2>(argc, argv);
    }
    break;
    case 3:
    {
      LabelGeometryMeasures<3>(argc, argv);
    }
    break;
    default:
      std::cout << "Unsupported dimension" << std::endl;
      return EXIT_FAILURE;
  }
  return EXIT_SUCCESS;
}
} // namespace ants
